1. Field of Invention
The invention describes a method utilizing actively illuminated imaging polarimetry for identifying a feature or features in an object through evaluation of how the object scrambles different incident fully polarized states into partially polarized states. Such a reduction in the degree of polarization of the light can be due to variations in the object below the resolution imaging system.
2. Description of the Background Art
An active polarization image of an object, usually in the form of a Mueller matrix image, contains pixel-by-pixel information on the intensity, the retardance, and the diattenuation (partial polarization) of the object at the spatial resolution of the camera or imaging system. When the proper measurements are performed and analyzed by the means of the present invention, information is obtained on the variations of the polarization properties within a pixel. Although imaging information regarding the spatial location of features within a pixel cannot be obtained below the resolution of the instrument, useful information is obtained on the variations of structures, the types of those structures, their orientations, and overall order when these structures modulate polarization spatially. The exploitation of these polarization parameters is useful in medical diagnostics (e.g., retinal disease diagnostics), crystal analysis, surveillance and many other applications.
The eye is subject to a number of pathologies and disease states which, if unchecked, lead to the impairment and possible loss of vision. Non-invasive optical methods are a preferential method for assessing the status of the eye because of ease of measurement. But the common eye diseases, including glaucoma, diabetic retinopathy, and Age-related Macular Degeneration (AMD) do not produce easily detectable changes during their early stages because many ocular structures are transparent and because others are obscured by opaque tissues. Improved optical ophthalmic diagnostics hold the promise of detecting ocular disease earlier and of quantifying the state of the disease more accurately so that its progression can be tracked and the effectiveness of treatments such as drug treatments can be evaluated.
One principal ophthalmic diagnostic is fundus imaging which takes a conventional image of the inside of the eye. To produce an image requires a contrast mechanism, in the case of the eye a distribution of chromophores. But the eye is nearly transparent from the cornea until the retinal pigment epithelium is encountered. In the back of the eye the dominant chromophores are hemoglobin in the retinal blood vessels and melanin in the retinal pigment epithelium (RPE). Thus there are not that many strong contrast mechanisms for imaging to exploit. Fundus imaging can detect neovascularization associated with diabetic retinopathy and the cupping of the optic nerve head associated with advanced glaucoma. But fundus imaging cannot detect AMD, early stage glaucoma, or the precursors to diabetic retinopathy.
Another successful ophthalmic imaging diagnostic is optical coherence tomography (OCT) which constructs a three dimensional image of retinal tissue by scanning a broad band point source while rapidly adjusting an interferometric delay line. OCT reveals details which fundus imaging cannot include the nerve fiber layer and structures deeper in the retina including the RPE, Bruchs membrane, and the choroid. Being interferometric, OCT images are speckle patterns making quantitative image analysis of small features very difficult. OCT does not measure the polarization of transparent birefringent structures or depolarization and thus misses much useful information on the retina.
To obtain additional information on the retina and other ocular structures polarization can be used to induce contrast in transparent tissues or to modify the contrast in tissues with chromophores. Polarization provides an additional set of contrast mechanisms through the three basic polarization properties, retardance, diattenuation, and depolarization. All light/matter polarization effects can be classified into these three groups.
The largest polarization effects in the eye are associated with retardance, the phase delay between polarization states which accumulates when light propagates through birefringent materials. In the eye large retardance is associated with the cornea and the nerve fiber layer and small but nonzero retardance with all other transparent tissues.
The present invention addresses methods of measuring other polarization effects and applying these methods to retinal imaging and other problems involving active polarimetric imaging. To prepare for this complex discussion of depolarization imaging, first the present ophthalmic polarization measuring technologies will be reviewed and their limitations considered, including the techniques of nerve fiber layer (NFL) retardance imaging and Mueller matrix imaging.
The NFL layer consists of retinal ganglion cell axons which are arranged into bundles of parallel fibers. This assembly of fibers is modeled as arrays of parallel, weakly reflecting nonabsorbing, dielectric cylinders embedded in a medium of slightly lower refractive index. The cylinder array model has been used to predict that the reflection of the NFL should be proportional to thickness, and that the backscattering reflection should be into a cone. The model has also been used to predict that the NFL possesses ‘linear form birefringence’ and behaves as a positive uniaxial crystal with optic axis parallel to the axis of the fibers. These predictions have been experimentally verified with some success; in particular, retardance magnitude has been shown to correlate well with NFL thickness. With more careful histology of retinal samples, it has been shown that a number of coaxial quasicylindrical structures of various diameters are present. As a result, the array model has been extended to include both thin cylinders (diameter ˜ 1/10 of a wavelength) and thick cylinders (diameter ˜1 wavelength). Considerable success has been achieved in using this model to explain experimentally obtained polarization measurements. Reported results indicate that the retina is a linear retarder, with retardance dependent on thickness; and that there is weak diattenuation and little depolarization on reflection. Reported retardance values are in the range of 2-4nm for NFL thickness of ˜15μm, corresponding to approximately 0.2nm/ μm.
Several obstacles exist to determining retinal polarization by direct measurement. Probably the most significant is that the polarization state of light entering the eye may be modified by other ocular structures, such as the cornea, crystalline lens, and vitreous. The vitreous and lens appear to have no significant polarization effect, but the cornea can have very substantial effect. Studies in humans suggest that the cornea is a linear retarder with slow axis pointing nasally downward, with the retardance magnitude increasing radially. The orientation and magnitude of the retardance has been shown to vary significantly between individuals, with double-pass retardances of 0-250 nm reported. Some researchers have attempted to avoid this problem by excising the anterior segment of the eye, and others have attempted to compensate for the corneal retardance by inserting a fixed or variable compensatory retarder in the measurement path.
A polarimeter for measurements of the NFL is commercially available from Laser Diagnostics, San Diego, Calif. and is called the GDx Retinal Polarimeter. The GDx is approved by the FDA for the diagnosis of glaucoma. The GDx is marketed as an instrument which assists in early glaucoma detection. The GDx is described in U.S. Pat. No. 6,704,106 to Anderson.
FIG. 4 in this patent is a diagram of the optical layout of the GDx and explains how the GDx acquires polarization information. A 780 nm laser diode is polarized and coupled into the main optical axis by a non-polarizing or nearly non-polarizing beam splitter. The polarization state is modulated by a rotating half wave linear retarder. A two axis scanner raster scans the beam. The light passes through an optional corneal polarization compensator which is only present is some models. A variable focus lens directs the light into the eye where the polarized light interacts with the tissues of the eye, in particular the cornea and retinal nerve fiber layer. The variable focus lens is adjusted to focus to a small spot of light on the retina. The raster scanned beam reflects and scatters from the back of the eye and the returned light constitutes an image. After each image the half wave linear retarder is rotated to change the generated polarization state and the analyzed polarization state. The light scattered from the eye passes back through the rotating half wave linear retarder and is divided at the non-polarizing beam splitter. The half of the beam reflected back toward the source is lost. The transmitted fraction enters a polarizing beam splitter where the s- and p-polarized components are split and detected by separate detectors. One detectors signal corresponds to the co-polarized signal, it is polarized parallel to the illuminating beam and is the portion of the light whose polarization was not changed. The other detector measures the cross-polarized signal, that component whose polarization state was changed. The GDx illuminates the eye with twenty linearly polarized beams each oriented 18° from the next generated by a rotating half wave linear retarder.
Thus, with the standard GDx, 20 images are acquired, ten co-polarized and ten cross-polarized. Pixel by pixel, from a sinusoidal fit to these intensity values, the linear retardance and the retardance orientation are determined which when combined constitute a retardance image of the eye. The GDx compensates for the retardance of the cornea to produce a retardance image of the retina from which the thickness of the nerve fiber layer is determined. A normal retardance distribution indicates NFL health. A thin NFL is an indication of glaucoma. A series of measurements can document the progressive thinning and estimate the rate of glaucoma progression.
The GDx as produced by Laser Diagnostics is an incomplete polarimeter. The GDx cannot measure a complete Mueller matrix, cannot measure diattenuation, and cannot measure any depolarization parameters as the present invention can. The GDx is designed to measure linear retardance, which correlates with the thickness and health of the NFL. The GDx does not measure the other polarization properties: diattenuation and depolarization. Thus, the GDx device does not provide a full set of depolarization images of the retina and, therefore, does not provide a complete set of information about the deeper layers of the retina, lesions associated with neovascularization and basal and basal laminar deposits.
A simplified version of FIG. 4 of U.S. Pat. No. 6,704,106 is shown in FIG. 1, where a laser emits light through a polarizer onto a non-polarizing beam splitter. One part of the light beam is reflected to the eye via a rotating half-wave retarder and an objective lens. Another part of the light beam is reflected toward a polarizing beam splitter. The polarizing beam splitter reflects co-polarized light to one photodetector and, optionally, cross-polarized light to another photodetector. Light reflected from the eye also passes through the two beam splitters to the photodetectors.
The aforementioned limitations are partially addressed by conventional Mueller matrix imaging of ocular tissues. However, before describing these conventional methods and systems for Mueller matrix imaging of ocular tissues, a background discussion of Mueller matrices is presented.
Depolarization can be calculated by comparing the input and output Degree of Polarization (DoP) of Stokes vectors. A more general method is via algorithms which operate on Mueller matrices. This is the subject of the present invention. These methods are more general because different polarization states are depolarized by different amounts, and these variations occur in many different ways, i.e. with different degrees of freedom. Before presenting the present inventions several quantities regarding Mueller matrices need to be given precise definitions. The present invention is a set of new algorithms to exploit these depolarization variations described by the Mueller matrix to infer additional information on the structure and subpixel order of the sample.
Frequent reference is made here to the description of polarization states by Stokes vectors, their location on the Poincaré sphere, and the properties of Mueller matrices, all defined and extensively discussed in 3.1.6 and 3.2.1 of Brossseau (1998) and many other standard works on optical polarization. All fully polarized Stokes vectors can be described by their location on the Surface of the Poincaré sphere by the equation
                              S          ⁡                      [                          θ              ,              ϕ                        ]                          =                              [                                                            1                                                                                                                        Cos                      ⁡                                              [                                                  2                          ⁢                          θ                                                ]                                                              ⁢                                          Cos                      ⁡                                              [                        ϕ                        ]                                                                                                                                                                                    Sin                      ⁡                                              [                                                  2                          ⁢                          θ                                                ]                                                              ⁢                                          Cos                      ⁡                                              [                        ϕ                        ]                                                                                                                                                              Sin                    ⁡                                          [                      ϕ                      ]                                                                                            ]                    =                                    [                                                                                          S                      0                                                                                                                                  S                      1                                                                                                                                  S                      3                                                                                                                                  S                      4                                                                                  ]                        .                                              (        1        )            where the parameter θis the orientation of polarization ellipse major axis and φ is the latitude on the Poincaré sphere.
The Degree of Polarization (DoP) of a Stokes vector is defined as
                              DoP          ⁡                      [            S            ]                          =                                                                              S                  1                  2                                +                                  S                  2                  2                                +                                  S                  3                  2                                                                    S              0                                .                                    (        2        )            
When the DoP of the exiting beam is less than the DoP of the incident beam, then depolarization occurred. In particular when the incident beam is completely polarized (DoP=1) and the exiting beam has a DoP=a, where 0≦a≦1, then the reduction in DoP, 1−a, indicates the depolarization that particular polarization state.
Two other Stokes vector parameters will be used later in the DoP maps. The orientation of the major axis of the polarization ellipse, θ, is given by
                    θ        =                              1            2                    ⁢          Arc          ⁢                                          ⁢                                    Tan              ⁡                              (                                                      S                    2                                    /                                      S                    1                                                  )                                      .                                              (        3        )            The Degree of Circular Polarization, DoCP, is
                              DoCP          ⁡                      [            S            ]                          =                                            S              3                                      S              0                                .                                    (        4        )            
A linear interaction of incident light with a sample has its polarization transformations properties are described by a Mueller matrix, M, that relates the incident Stokes vector, Sincident, with the exiting Stokes vector, SExiting, by the relation
                              S          Exiting                =                              M            ·                          S              Incident                                =                                                    [                                                                                                    M                        00                                                                                                            M                        01                                                                                                            M                        02                                                                                                            M                        03                                                                                                                                                M                        10                                                                                                            M                        11                                                                                                            M                        12                                                                                                            M                        13                                                                                                                                                M                        20                                                                                                            M                        21                                                                                                            M                        22                                                                                                            M                        23                                                                                                                                                M                        30                                                                                                            M                        31                                                                                                            M                        32                                                                                                            M                        33                                                                                            ]                            ·                              [                                                                                                    S                        0                                                                                                                                                S                        1                                                                                                                                                S                        2                                                                                                                                                S                        3                                                                                            ]                                      =                                          [                                                                                                    S                        0                        ′                                                                                                                                                S                        1                        ′                                                                                                                                                S                        2                        ′                                                                                                                                                S                        3                        ′                                                                                            ]                            .                                                          (        5        )            
Those skilled in the optics understand that nonlinear interactions of light with samples, such as occur at higher light intensity levels are described by similar but more complex equations where a different Mueller matrix is necessary for each generated wavelength.
The exiting DoP for each possible incident Stokes vector is given by the equation
                              DoP          ⁡                      (                          M              ,              S                        )                          =                                                                                                                        S                      1                      ′                                        ⁡                                          (                                              M                        ,                        S                                            )                                                        2                                +                                                                            S                      2                      ′                                        ⁡                                          (                                              M                        ,                        S                                            )                                                        2                                +                                                                            S                      3                      ′                                        ⁡                                          (                                              M                        ,                        S                                            )                                                        2                                                                                    S                0                ′                            ⁡                              (                                  M                  ,                  S                                )                                              .                                    (        6        )            This equation contains information on all the possible depolarizations which an optical element or polarimetry sample performs at the wavelength, angle of incidence, aperture, and other optical beam parameters where the Mueller matrix was measured, calculated, or simulated.
Before discussing the more complex depolarization properties which are the subject of the invention, several conventional polarization properties associated with a Mueller matrix shall be defined in the exact mathematical form which will be used in the following algorithms including the following: diattenuation, diattenuation vector, polarizance, polarizance vector, retardance, and retardance vector. These comprise the basic nondepolarizing properties of a polarization transformation.
Diattenuation
Diattenuation is the property of polarizers and partial polarizers whereby the transmission is a function of the incident polarization state. Diattenuation is entirely described by the first row of the Mueller matrix by elements m00, m01, m02, and m03. The diattenuation (also “diattenuation magnitude”), D, is
                              D          =                                                                                          m                    01                    2                                    +                                      m                    02                    2                                    +                                      m                    03                    2                                                                              m                00                                      =                                                            T                  max                                -                                  T                  min                                                                              T                  max                                +                                  T                  min                                                                    ,                            (        7        )            which describes the degree to which M is a partial polarizer, given by the maximum and minimum transmittance. Further these elements define two additional degrees of freedom on the Poincaré sphere, the diattenuation axis orientation and latitude as
                              (                                                    (                                  2                  ⁢                  θ                                )                            axis                        ,                          ϕ              axis                                )                =                              (                                          Arc                ⁢                                                                  ⁢                                  Cos                  [                                                            m                      01                                                                                                                m                          01                          2                                                +                                                  m                          02                          2                                                                                                      ]                                            ,                              Arc                ⁢                                                                  ⁢                                  Sin                  [                                                            m                      3                                                                                                                m                          01                          2                                                +                                                  m                          02                          2                                                +                                                  m                          03                          2                                                                                                      ]                                                      )                    .                                    (        8        )            
The diattenuation axis passes through the incident polarization states of maximum transmittance, Smax, also called the Diattenuation Vector,Smax=(√{square root over (m012+m022+m032)}, m01, m02, m03)  (9)and the incident polarization state of minimum transmittance, Smin,Smin=(√{square root over (m012+m022+m032)}, −m01, −m02, −m03)  (10)
Polarizance
Polarizance describes the ability of a polarization element to polarize unpolarized light. Polarizance, P, is entirely described by the first column of the Mueller matrix by elements m00, m10, m20, and m30, and is
                              P          =                                                                      m                  10                  2                                +                                  m                  20                  2                                +                                  m                  30                  2                                                                    m              00                                      ,                            (        11        )            which is the degree of polarization of the exiting light when unpolarized light is incident on Mueller matrix, M. The Polarizance Vector, PV, is this exiting state,PV=(√{square root over (m102+m202+m302)}, m10, m20, m30)  (12)
The Depolarization Index
A Mueller matrix which is not depolarizing is called a Non-depolarizing Mueller Matrix and satisfies the condition
                              4          ⁢                      M            00            2                          =                              ∑                          i              =              0                        3                    ⁢                                    ∑                              j                =                0                            3                        ⁢                                          M                                  i                  ,                  j                                2                            .                                                          (        13        )            
This constraint was manipulated by Gil and Bernabeu in 1985 to obtain the Depolarization Index, DI, the only single number metric in common usage for characterizing the depolarization of a Mueller matrix. The DI is defined as
                              DI          ⁡                      [            M            ]                          =                                                                              (                                                            ∑                                              i                        ,                                                  j                          =                          0                                                                    3                                        ⁢                                          M                                              i                        ,                        j                                            2                                                        )                                -                                  M                                      0                    ,                    0                                    2                                                                                    3                            ⁢                              M                                  0                  ,                  0                                                              .                                    (        14        )            The Depolarization Index” is in widespread use to characterize the “strength” or “magnitude” of the depolarization but we will see that this metric can produce misleading results. The DI can vary between zero and one. The DI equals one for non-depolarizing Mueller matrices and equals zero for the ideal depolarizer Mueller matrix, ID
                    ID        =                              [                                                            1                                                  0                                                  0                                                  0                                                                              0                                                  0                                                  0                                                  0                                                                              0                                                  0                                                  0                                                  0                                                                              0                                                  0                                                  0                                                  0                                                      ]                    .                                    (        15        )            the Mueller matrix which completely depolarizes all incident polarization states so that only unpolarized light exits. The numerator of the Depolarization Index is the Euclidean distance from the ideal depolarizer to the Mueller matrix in a sixteen dimensional space formed from the sixteen elements. The denominator is the radius of the hyperspherical surface of nondepolarizing Mueller matrices for that M00. So the DI is the fractional distance of a Mueller matrix along a line segment from ID to the hyperspherical surface for nondepolarizing Mueller matrices. The Depolarization Index is not a consistent metric of the amount of depolarization.
Problems with the Depolarization Index
Any matrix of the form
                              CD          =                      [                                                            1                                                  a                                                  b                                                  c                                                                              0                                                  0                                                  0                                                  0                                                                              0                                                  0                                                  0                                                  0                                                                              0                                                  0                                                  0                                                  0                                                      ]                          ,                                                            a                2                            +                              b                2                            +                              c                2                                              ≤          1.                                    (        16        )            is a complete depolarizer and only outputs unpolarized light with DoP=0 for any incident polarization state. But the Depolarization Index for these matrices varies over the range from zero to 0.577,
                              DI          =                                                                      a                  2                                +                                  b                  2                                +                                  c                  2                                                                                    3                            ⁢                              M                00                                                    ,                  0          ≤          DI          ≤                      1                          3                                ≈                      0.577            .                                              (        17        )            Therefore the Depolarization Index is not useful for describing the DoP output by Mueller matrices which can only output unpolarized light with DoP═0. DoP may have meaning in a sixteen-dimensional Mueller matrix space, but it has problems in describing the average DoP output by a highly depolarizing process with diattenuation, indicated by the parameters a, b, and c. Thus the present invention describes more useful and understandable metrics for describing depolarization than the Depolarization Index.
With the previous discussion of Mueller matrices as background, we may discuss U.S. Pat. No. 5,822,035 to Bille which describes a Mueller matrix ophthalmic polarimeter and its application to the measurement of corneal retardance. But despite measuring the Mueller matrix of the cornea and the ocular tissue Bille's invention does not describe how to determine the depolarization or diattenuation of ocular tissues or any other samples. Thus, Bille does not analyze determine variations of the depolarization with different incident polarization states of any material, let alone ocular tissue.
An example of retinal Mueller matrix imaging is described in the paper Confocal Scanning Laser Ophthahmoscopy of Mueller-matrix Polarimetry by Bueno and Campbell. They describe an optical system to perform Mueller matrix imaging measurement of the retina and provides an example of a retinal Mueller matrix image. They uses the Mueller matrix image to compare pixels in small regions to calculate the combination of illuminating polarization state and analyzer which provide the best and worst signal to noise ratio. They then apply the best combination to acquire images with improved image quality. But their method does not address how to characterize the depolarization or diattenuation of any material, let alone ocular tissue.
Depolarization of the retina is addressed in the very recent proceedings paper “Depth-Coded polarization imaging” by Elsner, Weber, Cheney, Smithwick, and Burns. They acquired images using the GDx by measuring forty intensities in the cross-polarized channel. These intensities have a roughly sinusoidal variation as a function of the rotating half wave retarder angle. The minimum intensity is taken as the most depolarized state and the intensity recorded pixel by pixel into an image. This measurement and algorithm can measure a mixture of depolarization, retardance, and diattenuation but as will be shown shortly is still several steps removed from a complete depolarization measurement. The method of Elsner et al. provides as close to a complete depolarization measurement as can be obtained with an off-the-shelf GDx. Further analysis of the method of Elsner et al. follows.
In the crossed polarized channel of the GDx the generated and analyzed polarizations are always orthogonal linear polarizations (90° apart). The generated signal comes from a linear polarized beam through a rotating half wave linear retarder producing a linear polarization which rotates at twice the speed of the retarder. The analyzed beam goes back through the same rotating half wave linear retarder then into a polarizer at 90° to the initial polarizer. The analyzed state is a rotating linear state 90° from the generated state.
Any depolarizer which is placed between crossed linear polarizers produces a depolarized (unpolarized) component which cannot be extinguished by the analyzing polarizer. With only the presumed linear retarder (the NFL) in the sample compartment a sinusoidal signal is produced which goes to zero intensity when the retarder is aligned with the polarization generator or polarization analyzer. The method of Elsner et al. exploits the fact that when the retarder also has a depolarization component the signal will not go to zero and the depolarization is sensed. Here, the minimum value in the sinusoidal signal is presumed to be due to depolarization.
One limitation of the method of Elsner et al. is that if the sample has retardance and diattenuation which are not aligned then the signal will not go to zero but the leakage is not due to depolarization. Another limitation of the method of Elsner et al. is that if the sample has elliptical diattenuation the signal will not go to zero and the depolarization measurement may not be accurate. A further limitation of the method of Elsner et al. is that if the sample has circular retardance or elliptical retardance the signal will not go to zero and the depolarization measurement may not be accurate. Thus, with the method of Elsner et al., is not possible to tell what fractions of depolarized light, light coupled by diattenuation, and light coupled by retardance are present in this polarization metric.
The method of Elsner et al. is a fine attempt to get close to the depolarization information in the sample but the ability to measure depolarization accurately is limited by the apparatus and algorithm. The present invention describes an apparatus and algorithm to accurately measure the depolarization, and since depolarization is a complex phenomena, to provide a plurality of metrics describing the samples depolarization.
Therefore what is desired, as discovered by the present inventors, is a polarization imaging method, apparatus, and computer program product which is not limited to measuring linear retardance but can also measure diattenuation and measure all nine degrees of freedom of depolarization. Such a device will be suited to detecting lesions in the retina and other parts of the eye, detecting the characteristics of retinal blood vessels, analyzing the order of the RPE cell packing, and many other imaging tasks. Such a device will also be suited to active imaging of the environment and remote sensing, as well as optical testing of optical systems, as well as other polarimetry operations.